MacHack 1994/MacHack 1994.toast/MacHack™ 1987-1994/MacHack™ '89/Papers '89/Horwat Papers/FactorLab

Text File | 1989-06-13 | 16.8 KB | 543 lines | [TEXT/CCL ]

;
;Prime-testing and factorization procedures
;
;Copyright 1987, 1989 Waldemar Horwat
;
;May 25, 1987.
;June 13, 1989.
;

;
;This program contains a few interesting routines.
;First there are several routines for adding, subtracting, multiplying, and
;exponentiating numbers in modular arithmetic, as well as implementing Euclid's
;algorithm
;
;The next major section is a set of primality testers.  The main entry point
;is the prime? function, which tests whether a number is prime using a number of
;tests.
;
;The following section consists of factoring routines.  Use the factor function
;to completely factor a number.
;
;Finally, a demo of the RSA algorithm is included.  This algorithm should not be
;used for encryption purposes without checking the RSA patent.
;

;This program is written in standard Common Lisp.  It works under Coral Common Lisp;
;it has not been tested on Pearl Lisp, which is a subset of Common Lisp.


;*************
;* Utilities *
;*************

(defun display* (&rest items)
  (dolist (item items) (prin1 item) (princ " ")))

(defun some? (f l)
  (declare (optimize speed))
  (cond ((endp l) nil)
        ((funcall f (car l)) t)
        (t (some? f (cdr l)))))


(proclaim '(inline square))
(defun square (x) (* x x))


;;;Return three values x, y, and g such that ax+by=g=gcd(a,b).
(defun euclid (a b)
  (defun euclid-sub (x1 x3 y1 y3)
    (if (zerop y3)
      (cond
       ((zerop x3) (values 0 0 0))
       ((< x3 0) (euclid-sub (- x1) (- x3) y1 y3))
       (t (values x1 (if (zerop b) 0 (/ (- x3 (* a x1)) b)) x3)))
      (multiple-value-bind (q r) (floor x3 y3)
        (euclid-sub y1 y3 (- x1 (* q y1)) r))))
  (euclid-sub 1 a 0 b))

;;; Compute the sum of the args modulo the modulus.
(defmacro m+ (modulus &rest args)
  `(mod (+ ,@args) ,modulus))

;;; Compute the difference or negation of the args modulo the modulus.
(defmacro m- (modulus &rest args)
  `(mod (- ,@args) ,modulus))

;;; Compute the product of the args modulo the modulus.
(defmacro m* (modulus &rest args)
  (cond ((null args) 1)
        ((endp (cdr args)) (car args))
        (t `(mod (* ,(car args) (m* ,modulus ,@(cdr args))) ,modulus))))

;;; Compute the square of n modulo the modulus.
(defun msquare (modulus n)
  (mod (* n n) modulus))

;;; Compute n^exp (mod modulus).  exp should be nonnegative.
(defun mexpt (modulus n exp)
  (cond ((zerop exp) 1)
        ((evenp exp)
         (msquare modulus (mexpt modulus n (floor exp 2))))
        (t (m* modulus n (square (mexpt modulus n (floor exp 2)))))))


(defun divisible? (n divisors)
  (some? #'(lambda (d) (zerop (rem n d))) divisors))



;*****************
;* Prime Testing *
;*****************


;Return a list containing the n smallest primes.
(defun n-smallest-primes (n)
  (defun next-n-primes (p n primes)
    (cond ((= n 0) primes)
          ((divisible? p primes) (next-n-primes (+ 2 p) n primes))
          (t (next-n-primes (+ 2 p) (1- n) (nconc primes (list p))))))
  (cond ((<= n 0) nil)
        ((= n 1) '(2))
        (t (next-n-primes 3 (- n 1) (list 2)))))

(defvar small-primes (n-smallest-primes 30))

;
;This is a quick primality test that checks that p contains no small factors.
;This test returns () if p is not a prime, t if it is, and maybe if it is
;likely a prime.
;

(defun prime0p (p)
  (prime0-check p small-primes))

(defun prime0-check (p primes)
  (cond ((< p 2) nil)
        ((endp primes) 'maybe)
        ((> (square (car primes)) p) t)
        ((zerop (rem p (car primes))) nil)
        (t (prime0-check p (cdr primes)))))



;
;This is a quick primality test that checks that x^p-1=1 (mod p), where p
;is the suspect prime and x is a small number.
;This test returns nil if p is not a prime, t if it is, and maybe if it is
;likely a prime.
;
;This test is no longer used because it might consistently fail on Carmichael numbers.
;

(defun prime1p (p)
  (or (= p 2)
      (= p 3)
      (= p 5)
      (and (> p 1)
           (= 1 (mexpt p 2 (1- p)))
           (= 1 (mexpt p 3 (1- p)))
           (= 1 (mexpt p 5 (1- p)))
           'maybe)))



;
;This primality test uses the Miller-Rabin test to check whether a number is
;prime.  This test may falsely assert that a number is prime with a probability
;no greater than 2^(-2*k), where k is the number of iterations.  If k=32, the
;probability of a wrong primality answer is no greater than 2^-64; actually, the
;probability is much smaller than this.
;k can be set by changing the *miller-iteration-count* global.
;

;Return nil if p is not a prime, and t if it is (may be wrong with probability
;as specified above).  Return maybe if this test has been turned off.

(defvar *miller-iteration-count* 32)

(defun prime2p (n)
  (cond
   ((<= *miller-iteration-count* 0) 'maybe)
   ((< n 2) nil)
   ((= n 2) t)
   ((evenp n) nil)
   (t (let ((minusone (- n 1)))
        (multiple-value-bind (s tt) (odd-shift minusone)
          (dotimes (i *miller-iteration-count* t)
            (let ((bb (mexpt n (1+ (random minusone)) tt)))
              (unless (or (= bb 1) (= bb minusone))
                (do ((b (msquare n bb) (msquare n b))
                     (j 1 (1+ j)))
                    ((= b minusone))
                  (when (or (= b 1) (> j s))
                    (return-from prime2p nil)))))))))))

;Return s and t such that n=t*2^s with t odd.
(defun odd-shift (n)
  (do ((tt n (/ tt 2))
       (s 0 (1+ s)))
      ((oddp tt) (values s tt))))

;Return t if the prime tests test1 and test2 are equivalent on p.
(defun equal-testp (p test1 test2)
  (eq (not (funcall test1 p)) (not (funcall test2 p))))

;Check that the prime tests test1 and test2 are equivalent over a range of
;numbers starting at n.
(defun equal-test-range (n test1 test2)
  (let ((t1 (funcall test1 n))
        (t2 (funcall test2 n)))
    (if t1 (display* n))
    (unless (eq (not t1) (not t2)) (cerror "Search will resume" "Fail on ~D" n))
    (equal-test-range (1+ n) test1 test2)))

(defvar prime-procedures (list #'prime0p #'prime2p))

;
;Test whether p is a prime.  Return true, false, or maybe.

(defun prime? (p)
  (prime-procedures? p prime-procedures))

(defun prime-procedures? (p procs)
  (if (endp procs) 'maybe
      (let ((answer (funcall (car procs) p)))
        (if (eq answer 'maybe)
          (prime-procedures? p (cdr procs))
          answer))))

;
;Find the smallest prime greater than or equal to n.
;

(defun next-prime (n)
  (cond ((<= n 2) 2)
        ((evenp n) (next-prime (1+ n)))
        ((prime? n) n)
        (t (next-prime (+ n 2)))))



;*****************
;* Factorization *
;*****************

;
;FactorList procedures.
;A FactorList is a partial factorization that consists of two lists:
;factors that are not necessarily primes and factors that are known to be
;prime.
;

;Create a FactorList with number n and no factors.
(proclaim '(inline new-primes-flist new-factors-flist add-prime-flist))
(defmacro new-primes-flist (&rest primes)
  (list 'list (cons 'list primes)))
(defmacro new-factors-flist (&rest numbers)
  (cons 'list (cons nil numbers)))

(defmacro make-flist (primes factors) `(cons ,primes ,factors))

(defmacro primes-flist (flist) `(car ,flist))
(defmacro factors-flist (flist) `(cdr ,flist))

(defun add-prime-flist (p flist)
  (make-flist (cons p (primes-flist flist))
              (factors-flist flist)))

(defun merge-lists (list1 list2)
  (cond ((endp list1) list2)
        ((endp list2) list1)
        ((<= (car list1) (car list2))
         (cons (car list1) (merge-lists (cdr list1) list2)))
        (t (cons (car list2) (merge-lists list1 (cdr list2))))))

(defun merge-flists (flist1 flist2)
  (make-flist (merge-lists (primes-flist flist1) (primes-flist flist2))
              (merge-lists (factors-flist flist1) (factors-flist flist2))))

;
;This factoring program attempts to factor n by checking whether it is a
;prime.  If so, it flags it as such.
;

(defun factor0 (n)
  (if (prime? n) (new-primes-flist n)
      (new-factors-flist n)))

;
;This is a quick factoring program that factors small primes out of p.
;This test returns a FactorList.
;

(defun factor1-primes (n primes)
  (if (endp primes) (new-factors-flist n)
      (let ((prime (car primes)))
        (cond ((> (square prime) n) (new-primes-flist n))
              ((zerop (rem n prime))
               (add-prime-flist prime (factor1-primes (floor n prime) primes)))
              (t (factor1-primes n (cdr primes)))))))

(defun factor1 (n)
  (if (= n 1) (new-factors-flist n)
      (factor1-primes n small-primes)))


;
;This is a Pollard rho factoring program.
;

(defvar *visible-rho-factor* t)

(defun rho-factor1 (n) (rho-factor n 1))
(defun rho-factor3 (n) (rho-factor n 3))

(defun rho-factor (n coeff)
  (if (< n 4) (new-factors-flist n)
      (rho-factor-sub n coeff 1 1 1 1 2 256)))

(defun rho-factor-sub (n coeff x1 x2 c1 c2 climit step)
  (defun iterate (x1 xref count)
    (do ((product 1 (rem (* product (abs (- x xref)) product) n))
         (x x1)
         (i count (1- i)))
        ((zerop i) (values x (gcd product n)))
      (setq x (rem (+ (square x) coeff) n))))
  (if *visible-rho-factor* (format t "[~D] " c1))
  (let ((count (min (- climit c1) step)))
    (multiple-value-bind (x g) (iterate x1 x2 count)
      (cond ((= 1 g)
             (if (= climit (+ c1 count))
               (rho-factor-sub n coeff x x climit climit (+ climit climit) step)
               (rho-factor-sub n coeff x x2 (+ c1 count) c2 climit step)))
            ((= n g)
             (if (= step 1)
               (new-factors-flist n)
               (rho-factor-sub n coeff x1 x2 c1 c2 climit (max 1 (floor step 16)))))
            (t (new-factors-flist g (floor n g)))))))


(defvar factor-procedures (list #'factor1 #'factor0 #'rho-factor1 #'rho-factor3))

;;; Factor n and return a list of factors.
(defun factor (n)
  (let ((l (factor-flist n)))
    (if (endp (cdr l)) (car l)
        (cerror "" "Unable to factor ~D: ~S" n l))))

;;; Factor n and return a factorization of n consisting of either lone factors
;;; or sublists of the form (p e) to indicate that p^e is a factor of n.  The
;;; factors are in increasing order and are not repeated.

(defun factorlist->ilist (factorlist)
  (if (null factorlist) ()
      (do
        ((first (car factorlist))
         (fl (cdr factorlist) (cdr fl))
         (count 1 (1+ count)))
        ((or (endp fl) (/= first (car fl)))
         (if (= count 1) (cons first (factorlist->ilist fl))
             (acons first count (factorlist->ilist fl)))))))

(defun factor-ilist (n)
  (factorlist->ilist (factor n)))

;
;Factor n and return a FactorList.
;

(defun factor-flist (n)
  (procs-factor n factor-procedures))

(defvar *visible-factor* t)

(defun procs-factor (n procs)
  (if (endp procs) (new-factors-flist n)
      (progn
        (if *visible-factor*
          (format t "Factoring ~D using ~S... " n (car procs)))
        (let ((factoring (funcall (car procs) n)))
          (if *visible-factor*
            (format t "~S~%" factoring))
          (cond ((null (factors-flist factoring)) factoring)
                ((= (car (factors-flist factoring)) n)
                 (procs-factor n (cdr procs)))
                (t (merge-flists
                    (make-flist (primes-flist factoring) nil)
                    (reduce #'merge-flists
                            (mapcar #'factor-flist (factors-flist factoring))))))))))

(defun fast-factor-sub (n k primes)
  (let ((prime-limit (square (next-prime (car (last primes))))))
    (defun check (n)
      (defun factor-power (prime)
        (if (zerop (rem n prime))
          (progn
            (setq n (floor n prime))
            (1+ (factor-power prime)))
          0))
      (let ((powers (mapcar #'factor-power primes)))
        (if (< n prime-limit)
          (list powers n)
          nil)))
    (let* ((r (floor (sqrt n)))
           (2r (* 2 r))
           (new-v)
           (new-a))
      (format t "Factor: n=~D k=~D r=~D primes=~S~%" n k r primes)
      (do ((old-u 2r u)
           (u 2r (- 2r (mod u new-v)))
           (old-v (- (* n k) (square r)) v)
           (v 1 new-v)
           (a 0 new-a)
           (old-p 1 p)
           (p r (mod (+ old-p (* new-a p)) n))
           (s 1 (- s))
           (count 0 (1+ count)))
          (nil)
        (format t "Iteration ~D: u=~D v=~D a=~D p=~D s=~D " count u v a old-p s)
        (prin1 (check v))
        (terpri)
        (setq new-v (+ old-v (* a (- old-u u))))
        (setq new-a (floor u new-v))))))

(defun phi (n) (phi-ilist (factor-ilist n)))
(defun phi-ilist (ilist)
  (cond ((null ilist) 1)
        ((numberp (car ilist)) (* (1- (car ilist)) (phi-ilist (cdr ilist))))
        (t (* (1- (caar ilist)) (expt (caar ilist) (1- (cdar ilist)))
              (phi-ilist (cdr ilist))))))

(defun jacobi (num den)
  (defun internal-jacobi (num den factor)
    (cond
     ((>= num den) (internal-jacobi (mod num den) den factor))
     ((= num 1) factor)
     ((evenp num) (internal-jacobi (/ num 2) den
                                   (if (eq (logbitp 1 den) (logbitp 2 den))
                                     factor
                                     (- factor))))
     ((evenp den) (internal-jacobi num (/ den 2) factor))
     (t (internal-jacobi den num
                         (if (and (logbitp 1 num) (logbitp 1 den))
                           (- factor)
                           factor)))))
  (cond
   ((/= (gcd num den) 1) 0)
   ((< num 0) (internal-jacobi (mod num den) den 1))
   (t (internal-jacobi num den 1))))

(defun jacobi-check (num den)
  (eq (mod (jacobi num den) den) (mexpt den num (/ (1- den) 2))))

(defun generator? (g p &optional p-ilist)
  (defun generator-check (g p p-ilist)
    (cond
     ((null p-ilist) t)
     ((= (mexpt p g (/ (1- p) (if (numberp (car p-ilist)) (car p-ilist) (caar p-ilist)))) 1)
      nil)
     (t (generator-check g p (cdr p-ilist)))))
  (cond
   ((null p-ilist) (generator? g p (factor-ilist (1- p))))
   ((/= (gcd g p) 1) nil)
   (t (generator-check g p p-ilist))))

(defun next-generator (g p &optional p-ilist)
  (cond
   ((null p-ilist) (next-generator g p (factor-ilist (1- p))))
   ((generator? g p p-ilist) g)
   (t (next-generator (1+ g) p p-ilist))))

(defun next-prime-with-gen (g p)
  (let ((prime (next-prime p)))
    (if (generator? g prime) prime
        (next-prime-with-gen g (1+ prime)))))

;(factor 3463544557620305157)
;(factor 80209029642225101201) ==> (8872758161 9039920641) [150784 iterations]


(defvar *visible-rho-ind* t)

(defun rho-ind (a g p)
  (let ((low-third (floor p 3))
        (high-third (floor (* 2 p) 3))
        (phi (1- p)))
    (defun rho (xlist)
      (let ((x (car xlist)))
        (cond
         ((<= x low-third)
          (list (m* p a x) (mod (1+ (cadr xlist)) phi) (caddr xlist)))
         ((<= x high-third)
          (list (msquare p x) (m* phi 2 (cadr xlist)) (m* phi 2 (caddr xlist))))
         (t
          (list (m* p g x) (cadr xlist) (mod (1+ (caddr xlist)) phi))))))
    (do
      ((x (rho '(1 0 0)) (rho x))
       (xn 1 (1+ xn))
       (y (rho (rho '(1 0 0))) (rho (rho y)))
       (yn 2 (+ 2 yn)))
      ((progn
         (if *visible-rho-ind*
           (format t "x~D = ~D = a^~D*g^~D, x~D = ~D = a^~D*g^~D~%"
                   xn (car x) (cadr x) (caddr x)
                   yn (car y) (cadr y) (caddr y)))
         (= (car x) (car y)))
       (let ((aexp (- (cadr x) (cadr y)))
             (gexp (- (caddr y) (caddr x))))
         (format t "a^~D = g^~D~%" aexp gexp)
         (multiple-value-bind (b c aexp) (euclid aexp phi)
           (let ((gexp (m* phi gexp b)))
             (format t "a^~D = g^~D~%" aexp gexp))))))))


;*******
;* RSA *
;*******

;Find a prime number using start as a true random number seed.
;The prime returned is congruent to 2 modulo 3.

(defun find-rsa-prime (start)
  (do ((p (+ 5 (* start 6)) (+ p 6)))
      ((prime? p) p)))

;Generate an RSA encoding key-decoding key pair.
;The two inputs should be random numbers from which the prime number search
;will begin.  Start1 and start2 should have approximately half as many digits
;as E and D.  Start1 and start2 should be different; the ratio of start1 and
;start2 should not be close to an exact fraction.  In a real application a few
;further precautions should be taken, as discussed in Knuth's volume 2.
;Return two values:  the modulus E and the decrypting key D.

(defun make-rsa (start1 start2)
  (let ((p (find-rsa-prime start1))
        (q (find-rsa-prime start2)))
    (let ((e (* p q)))
      (multiple-value-bind (x y g) (euclid 3 (* (1- p) (1- q)))
        (unless (= g 1)
          (error "3 and (p-1)(q-1) should be relatively prime."))
        (when (< x 0) ;Make x positive if necessary.
          (setq x (+ x (* (1- p) (1- q)))))
        (values e x)))))


;Encrypt the plaintext using the encrypting key E.

(defun encrypt-rsa (E plaintext)
  (when (< (square plaintext) E)
    (format t "Warning: plaintext is too short.~%"))
  (when (>= plaintext E)
    (format t "Warning: plaintext is too long; information was lost.~%"))
  (mexpt E plaintext 3))

;Decrypt the cyphertext using the decrypting key D.

;The encrypting key E is also needed.
(defun decrypt-rsa (E D cyphertext)
  (mexpt E cyphertext D))